
library(apsrtable)
library(sandwich)
library(xtable)
library(plyr)

load('pooled.frame.RData')

fit.w.robust <- function(model, cluster.var, dta){	
  robust.se <- function(model, cluster){
    require(sandwich)
    require(lmtest)
    M <- length(unique(cluster))
    N <- length(cluster)
    K <- model$rank
    dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
    uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
    rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
    rcse.se <- coeftest(model, rcse.cov)
    return(list(rcse.cov, rcse.se))
  }
  
  clustervar <- mapply(paste,cluster.var,dta[!(1:dim(dta)[1] %in% model$na.action),cluster.var],sep="")
  vcov <- robust.se(model, clustervar)
  model$se <- sqrt(diag(vcov[[1]]))
  model$coeftest <- vcov[[2]]
  return(model)
}

###
#Correct / Incorrect as Focus
###

pooled.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$category=="covid" & !is.na(pooled.frame$pid.binary)),])
pooled.model.covid.correct <- fit.w.robust(pooled.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$category=="covid" & !is.na(pooled.frame$pid.binary)),])

pooled.model.other.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$category!="covid" & !is.na(pooled.frame$pid.binary)),])
pooled.model.other.correct <- fit.w.robust(pooled.model.other.correct,'caseid',pooled.frame[which(pooled.frame$category!="covid" & !is.na(pooled.frame$pid.binary)),])

#Table 3
apsrtable(pooled.model.covid.correct,caption='Effect of Incentives on Pr(Correct Answer)',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('COVID-19 Information'),notes=list(stars.note,'Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="incentive|Intercept",coefnames, invert=TRUE)))

#Table D3
apsrtable(pooled.model.other.correct,caption='Effect of Incentives on Pr(Correct Answer)',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('COVID-19 Information'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="incentive|Intercept",coefnames, invert=TRUE)))

mortality.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$topic=="mortality" & !is.na(pooled.frame$pid.binary)),])
mortality.model.covid.correct <- fit.w.robust(mortality.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$topic=="mortality" & !is.na(pooled.frame$pid.binary)),])

international.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$topic=="international" & !is.na(pooled.frame$pid.binary)),])
international.model.covid.correct <- fit.w.robust(international.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$topic=="international" & !is.na(pooled.frame$pid.binary)),])

age.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$topic=="age" & !is.na(pooled.frame$pid.binary)),])
age.model.covid.correct <- fit.w.robust(age.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$topic=="age" & !is.na(pooled.frame$pid.binary)),])

source.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$topic=="source" & !is.na(pooled.frame$pid.binary)),])
source.model.covid.correct <- fit.w.robust(source.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$topic=="source" & !is.na(pooled.frame$pid.binary)),])

treatment.model.covid.correct <- lm(correct.answer ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$topic=="treatments" & !is.na(pooled.frame$pid.binary)),])
treatment.model.covid.correct <- fit.w.robust(treatment.model.covid.correct,'caseid',pooled.frame[which(pooled.frame$topic=="treatments" & !is.na(pooled.frame$pid.binary)),])

#Table C2
apsrtable(mortality.model.covid.correct,international.model.covid.correct,age.model.covid.correct,source.model.covid.correct,treatment.model.covid.correct,caption='Effect of Incentives on Pr(Correct Answer)',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('Mortality',"International","Risk by Age","Source","Treatment"),notes=list('Robust standard errors, clustered by Respondent, in parentheses'),omitcoef=expression(grep(pattern="incentive|Intercept",coefnames, invert=TRUE)))

#With Certainty Measure
pooled.model.covid.certainty <- lm(with.certainty ~ I(incentive %in% c('low')) + I(incentive %in% c('high')),data=pooled.frame[which(pooled.frame$category=="covid" & !is.na(pooled.frame$pid.binary)),])
pooled.model.covid.certainty <- fit.w.robust(pooled.model.covid.certainty,'caseid',pooled.frame[which(pooled.frame$category=="covid" & !is.na(pooled.frame$pid.binary)),])

#Table C5
apsrtable(pooled.model.covid.certainty,caption='Effect of Incentives on Pr(Correct Answer)',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('COVID-19 Information'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="incentive|Intercept",coefnames, invert=TRUE)))

interaction.model.correct <- lm(correct.answer ~ I(incentive %in% c('low'))*I(category=='covid') +  I(incentive %in% c('high'))*I(category=='covid'),data=pooled.frame)
interaction.model.correct <- fit.w.robust(interaction.model.correct,'caseid',pooled.frame)

#Table D6
apsrtable(interaction.model.correct,coef.names=c("(Intercept)","Low Incentive","COVID-19 Item","High Incentive","Low Incentive*COVID-19 Item","High Incentive*COVID-19 Item"))

###
#Partisan Divide as focus
###

pooled.model.covid.full <- lm(partisan.answer ~ I(pid.binary=='Democrat')*I(incentive %in% c('low')) + I(pid.binary=='Democrat')*I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category=="covid"),])
pooled.model.covid.full <- fit.w.robust(pooled.model.covid.full,'caseid',pooled.frame[which(pooled.frame$category=="covid"),])

pooled.model.other.full <- lm(partisan.answer ~ I(pid.binary=='Democrat')*I(incentive %in% c('low')) + I(pid.binary=='Democrat')*I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category!="covid"),])
pooled.model.other.full <- fit.w.robust(pooled.model.other.full,'caseid',pooled.frame[which(pooled.frame$category!="covid"),])

#Table C6
apsrtable(pooled.model.covid.full,caption='Effect of Incentives on Partisan Information Divide',coef.names=c("(Intercept)","Democrat","Low Incentive","High Incentive","Democrat*Low Incentive","Democrat*High Incentive"),model.names=c('COVID-19 Information'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="pid|incentive|Intercept",coefnames, invert=TRUE)))

#Table D5
apsrtable(pooled.model.other.full,caption='Effect of Incentives on Partisan Information Divide',coef.names=c("(Intercept)","Democrat","Low Incentive","High Incentive","Democrat*Low Incentive","Democrat*High Incentive"),model.names=c('Other Political Information'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="pid|incentive|Intercept",coefnames, invert=TRUE)))

#Baseline Correct/Incorrect
#Covid
descriptive.frame <- subset(pooled.frame,pooled.frame$incentive=="none" & !is.na(pooled.frame$pid.binary) & pooled.frame$category=='covid')

source.frame <- lm(correct.answer ~ pid.binary,data=descriptive.frame,subset=which(topic=="source"))
source.frame$coeftest <- coeftest(source.frame,vcov=vcovHC(source.frame,type="HC1"))

age.frame <- lm(correct.answer ~ pid.binary,data=descriptive.frame,subset=which(topic=="age"))
age.frame$coeftest <- coeftest(age.frame,vcov=vcovHC(age.frame,type="HC1"))

mortality.frame <- lm(correct.answer ~ pid.binary,data=descriptive.frame,subset=which(topic=="mortality"))
mortality.frame$coeftest <- coeftest(mortality.frame,vcov=vcovHC(mortality.frame,type="HC1"))

international.frame <- lm(correct.answer ~ pid.binary,data=descriptive.frame,subset=which(topic=="international"))
international.frame$coeftest <- coeftest(international.frame,vcov=vcovHC(international.frame,type="HC1"))

treatments.frame <- lm(correct.answer ~ pid.binary,data=descriptive.frame,subset=which(topic=="treatments"))
treatments.frame$coeftest <- coeftest(treatments.frame,vcov=vcovHC(treatments.frame,type="HC1"))

marginal.effects <- function(model){
  me.out <- model$coeftest
  dem.point <- me.out[1,1]
  dem.se <- me.out[1,2]
  
  diff.point <- me.out[2,1]
  diff.se <- me.out[2,2]
  
  rep.point <- dem.point + diff.point
  rep.se <- sqrt(vcov(model)[1,1] + vcov(model)[2,2] - 2*vcov(model)[2,1])

  point.vec <- c(dem.point,rep.point,diff.point)
  se.vec <- c(dem.se,rep.se,diff.se)
  
  point.diff <- cbind.data.frame(point.vec,se.vec)
  point.diff$lower <- point.diff[,1] - 2*point.diff[,2]
  point.diff$upper <- point.diff[,1] + 2*point.diff[,2]
  output <- point.diff
  return(output)
}

#Quantities for Table 1
marginal.effects(mortality.frame)
marginal.effects(source.frame)
marginal.effects(treatments.frame)
marginal.effects(international.frame)
marginal.effects(age.frame)

#Information Items on COVID with Certainty
answers.out.certainty <- ddply(descriptive.frame,.(pid.binary,topic),summarise,certainty=mean(with.certainty))
dem.answers <- answers.out.certainty[which(answers.out.certainty$pid.binary=='Democrat'),]
dem.answers$pid.binary <- NULL
names(dem.answers) <- c('topic','dem.correct')
rep.answers <- answers.out.certainty[which(answers.out.certainty$pid.binary=='Republican'),]
rep.answers$pid.binary <- NULL
names(rep.answers) <- c('topic','rep.correct')
party.divide.table <- merge(dem.answers,rep.answers,by=c('topic'),all.x=TRUE)
party.divide.table$partisan.difference <- abs(party.divide.table$dem.correct - party.divide.table$rep.correct)
party.divide.table <- party.divide.table[order(party.divide.table$partisan.difference,decreasing=TRUE),]
names(party.divide.table) <- c('topic','Share Correct (Dem)','Share Correct (Rep)','Partisan Divide')
rownames(party.divide.table) <- NULL
#Table C4
xtable(party.divide.table)

#Information in Table C3
table(pooled.frame$with.certainty[which(pooled.frame$category=="covid" & pooled.frame$incentive=='none' & pooled.frame$correct.answer==1 & !is.na(pooled.frame$pid.binary))]) #Higher values = More confidence among this group
table(pooled.frame$with.certainty[which(pooled.frame$category=="covid" & pooled.frame$incentive=='none' & pooled.frame$correct.answer==0 & !is.na(pooled.frame$pid.binary))]) #Higher values = Less Confidence among this group

#Information Items on Non-COVID Items
descriptive.frame.other <- subset(pooled.frame,pooled.frame$incentive=="none" & !is.na(pooled.frame$pid.binary) & pooled.frame$category!='covid')
answers.out.other <- ddply(descriptive.frame.other,.(pid.binary,topic),summarise,correct=mean(correct.answer))
dem.answers.other <- answers.out.other[which(answers.out.other$pid.binary=='Democrat'),]
dem.answers.other$pid.binary <- NULL
names(dem.answers.other) <- c('topic','dem.correct')
rep.answers.other <- answers.out.other[which(answers.out.other$pid.binary=='Republican'),]
rep.answers.other$pid.binary <- NULL
names(rep.answers.other) <- c('topic','rep.correct')
party.divide.table.other <- merge(dem.answers.other,rep.answers.other,by=c('topic'),all.x=TRUE)
party.divide.table.other$partisan.difference <- abs(party.divide.table.other$dem.correct - party.divide.table.other$rep.correct)
party.divide.table.other <- party.divide.table.other[order(party.divide.table.other$partisan.difference,decreasing=TRUE),]
names(party.divide.table.other) <- c('topic','Share Correct (Dem)','Share Correct (Rep)','Partisan Divide')
rownames(party.divide.table.other) <- NULL
#Table D1
xtable(party.divide.table.other)

####
##News Choice
####

pooled.frame$expert.choice <- ifelse(pooled.frame$news.choice=='Expert', 1, 0)
pooled.frame$mainstream.choice <- ifelse(pooled.frame$news.choice=='Mainstream',1,0)

pooled.frame$rightmedia.choice <- 0
pooled.frame$rightmedia.choice[which(pooled.frame$pid.binary=='Democrat' & pooled.frame$news.choice=='Outpartisan')] <- 1
pooled.frame$rightmedia.choice[which(pooled.frame$pid.binary=='Republican' & pooled.frame$news.choice=='Copartisan')] <- 1

pooled.frame$leftmedia.choice <- 0
pooled.frame$leftmedia.choice[which(pooled.frame$pid.binary=='Democrat' & pooled.frame$news.choice=='Copartisan')] <- 1
pooled.frame$leftmedia.choice[which(pooled.frame$pid.binary=='Republican' & pooled.frame$news.choice=='Outpartisan')] <- 1

newschoice.out <- ddply(pooled.frame[which(pooled.frame$incentive=='none' & pooled.frame$category=="covid"),],.(pid.binary),summarise,expert.choice=mean(expert.choice),mainstream.choice=mean(mainstream.choice),left.choice=mean(leftmedia.choice),right.choice=mean(rightmedia.choice))
newschoice.out <- as.data.frame(t(newschoice.out[c(1,2),c(2,3,4,5)]))
newschoice.out$divide <- abs(newschoice.out[,1] - newschoice.out[,2])
newschoice.out <- newschoice.out[order(newschoice.out$divide,decreasing=TRUE),]
names(newschoice.out) <- c('Share Selected (Dem)','Share Selected (Rep)','Partisan Divide')
#Table C1
xtable(newschoice.out)

newschoice.plot <- newschoice.out[,c(1,2)]
newschoice.plot <- t(as.matrix(newschoice.plot))

newschoice.frame <- subset(pooled.frame,!is.na(pooled.frame$pid.binary) & pooled.frame$category=="covid")

expert.choice <- lm(expert.choice ~ pid.binary,data=newschoice.frame)
expert.choice  <- fit.w.robust(expert.choice,'caseid',newschoice.frame)

mainstream.choice <- lm(mainstream.choice ~ pid.binary,data=newschoice.frame)
mainstream.choice  <- fit.w.robust(mainstream.choice,'caseid',newschoice.frame)

right.choice <- lm(rightmedia.choice ~ pid.binary,data=newschoice.frame)
right.choice  <- fit.w.robust(right.choice,'caseid',newschoice.frame)

left.choice <- lm(leftmedia.choice ~ pid.binary,data=newschoice.frame)
left.choice  <- fit.w.robust(left.choice,'caseid',newschoice.frame)

newsmarginal.effects <- function(model){
  me.out <- model$coeftest
  dem.point <- me.out[1,1]
  dem.se <- me.out[1,2]
  
  diff.point <- me.out[2,1]
  diff.se <- me.out[2,2]
  
  rep.point <- dem.point + diff.point
  rep.se <- sqrt(vcov(model)[1,1] + vcov(model)[2,2] - 2*vcov(model)[2,1])
  
  point.vec <- c(dem.point,rep.point)
  se.vec <- c(dem.se,rep.se)
  
  point.diff <- cbind.data.frame(point.vec,se.vec)
  point.diff$lower <- point.diff[,1] - 2*point.diff[,2]
  point.diff$upper <- point.diff[,1] + 2*point.diff[,2]
  output <- point.diff
  return(output)
}

right.barplot <- newsmarginal.effects(right.choice)
mainstream.barplot <- newsmarginal.effects(mainstream.choice)
expert.barplot <- newsmarginal.effects(expert.choice)
left.barplot <- newsmarginal.effects(left.choice)

newschoice.plot <- cbind.data.frame(right.barplot[,1],mainstream.barplot[,1],expert.barplot[,1],left.barplot[,1])
newschoice.plot <- as.matrix(newschoice.plot)

newschoice.plot.se <- cbind.data.frame(right.barplot[,2],mainstream.barplot[,2],expert.barplot[,2],left.barplot[,2])
newschoice.plot.se <- as.matrix(newschoice.plot.se)

error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
  arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}

#Figure 1
pdf(height=5,width=6,file="figure1.pdf")
par(mar=c(4.1,4,.5,.5))
mainplot <- barplot(newschoice.plot,beside=TRUE,ylim=c(0,.55),col=c('gray40','gray90'),las=1,names.arg=c("Right-Leaning","Mainstream","Health Expert","Left-Leaning"),ylab="Share Selecting Source",legend=c("Democrats","Republicans"),xlab="Media Source Type")
error.bar(mainplot,newschoice.plot,newschoice.plot.se)
box()
dev.off()

#Other Choices

newschoice.out.other <- ddply(pooled.frame[which(pooled.frame$incentive=='none' & pooled.frame$category!="covid"),],.(pid.binary),summarise,expert.choice=mean(expert.choice),mainstream.choice=mean(mainstream.choice),left.choice=mean(leftmedia.choice),right.choice=mean(rightmedia.choice))
newschoice.out.other <- as.data.frame(t(newschoice.out.other[c(1,2),c(2,3,4,5)]))
newschoice.out.other$divide <- abs(newschoice.out.other[,1] - newschoice.out.other[,2])
newschoice.out.other <- newschoice.out.other[order(newschoice.out.other$divide,decreasing=TRUE),]
names(newschoice.out.other) <- c('Share Selected (Dem)','Share Selected (Rep)','Partisan Divide')

#Table D2
xtable(newschoice.out.other)

###
#News Choice Pooled
###

#News Sources on COVID Information
pooled.model.covid.choice.copartisan <- lm(I(news.choice=='Copartisan') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category=="covid"),])
pooled.model.covid.choice.copartisan <- fit.w.robust(pooled.model.covid.choice.copartisan,'caseid',pooled.frame[which(pooled.frame$category=="covid"),])

pooled.model.covid.choice.expert <- lm(I(news.choice=='Expert') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category=="covid"),])
pooled.model.covid.choice.expert <- fit.w.robust(pooled.model.covid.choice.expert,'caseid',pooled.frame[which(pooled.frame$category=="covid"),])

pooled.model.covid.choice.mainstream <- lm(I(news.choice=='Mainstream') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category=="covid"),])
pooled.model.covid.choice.mainstream <- fit.w.robust(pooled.model.covid.choice.mainstream,'caseid',pooled.frame[which(pooled.frame$category=="covid"),])

pooled.model.covid.choice.outpartisan <- lm(I(news.choice=='Outpartisan') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category=="covid"),])
pooled.model.covid.choice.outpartisan <- fit.w.robust(pooled.model.covid.choice.outpartisan,'caseid',pooled.frame[which(pooled.frame$category=="covid"),])

#Table 2
apsrtable(pooled.model.covid.choice.copartisan,pooled.model.covid.choice.expert,pooled.model.covid.choice.mainstream,pooled.model.covid.choice.outpartisan,caption='Effect of Incentives on Information Source Selection',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('Copartisan News','Expert','Mainstream News','Out-Partisan News'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="pid|incentive|Intercept",coefnames, invert=TRUE)))

#News Sources on Other Information
pooled.model.other.choice.copartisan <- lm(I(news.choice=='Copartisan') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category!="covid"),])
pooled.model.other.choice.copartisan <- fit.w.robust(pooled.model.other.choice.copartisan,'caseid',pooled.frame[which(pooled.frame$category!="covid"),])

pooled.model.other.choice.expert <- lm(I(news.choice=='Expert') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category!="covid"),])
pooled.model.other.choice.expert <- fit.w.robust(pooled.model.other.choice.expert,'caseid',pooled.frame[which(pooled.frame$category!="covid"),])

pooled.model.other.choice.mainstream <- lm(I(news.choice=='Mainstream') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category!="covid"),])
pooled.model.other.choice.mainstream <- fit.w.robust(pooled.model.other.choice.mainstream,'caseid',pooled.frame[which(pooled.frame$category!="covid"),])

pooled.model.other.choice.outpartisan <- lm(I(news.choice=='Outpartisan') ~ I(incentive %in% c('low')) + I(incentive %in% c('high')) + factor(topic),data=pooled.frame[which(pooled.frame$category!="covid"),])
pooled.model.other.choice.outpartisan <- fit.w.robust(pooled.model.other.choice.outpartisan,'caseid',pooled.frame[which(pooled.frame$category!="covid"),])

#Table D4
apsrtable(pooled.model.other.choice.copartisan,pooled.model.other.choice.expert,pooled.model.other.choice.mainstream,pooled.model.other.choice.outpartisan,caption='Effect of Incentives on Information Source Selection',coef.names=c("(Intercept)","Low Incentive","High Incentive"),model.names=c('Copartisan News','Expert','Mainstream News','Out-Partisan News'),notes=list('Robust standard errors, clustered by Respondent, in parentheses','Models includes Topic Fixed Effects'),omitcoef=expression(grep(pattern="pid|incentive|Intercept",coefnames, invert=TRUE)))